devtools::install_github("kosukeimai/redist", ref = "dev")

library(tidyverse)
library(maptools)
library(igraph)
library(parallel)
library(redist)
library(spdep)
library(gridExtra)
library(sf)

ipw <- function(x, pop = NULL, com = NULL){
    if(!is.null(com)){
        indcom <- which(x$compact_ssd <= com)
    }else{
        indcom <- 1:ncol(x$partitions)
    }
    if(!is.null(pop)){
        indpop <- which(x$distance_parity <= pop)
    }else{
        indpop <- 1:ncol(x$partitions)
    }
    indbeta <- which(x$beta_sequence == 1)
    inds <- Reduce(intersect, list(indpop, indcom, indbeta))
    psi <- x[["energy_psi"]][inds]
    weights <- 1/exp(-1 * psi)
    inds <- sample(inds, length(inds), replace = TRUE, prob = weights)
    x <- x$partitions[,inds]
    return(x)
}
ben_theme <- function(){
    theme_classic() +
        theme(panel.background = element_blank(),
              panel.grid.major = element_blank(),
              panel.grid.minor = element_blank(),
              axis.line = element_line(colour = "black"),
              panel.border = element_rect(colour = "black", fill = NA, size = 1),
              strip.background = element_blank(),
              legend.position = "bottom", legend.title = element_blank(),
              plot.title = element_text(hjust = 0.5),
              plot.subtitle = element_text(hjust = .5))
}

fl_map <- readShapePoly("../data/fl/FL.shp")

## Load nodes
mn <- 89 
nodes <- strsplit(readLines("../data/largemap_enum_all/map_sub_nodes_89.dat")," ")[[1]]

## First, subset down to right map
fl_sub <- fl_map[which(rownames(fl_map@data) %in% nodes),]

## Then, reorder nodes to be correct order to match - otherwise, won't
## line up properly with original map
fl_sub <- fl_sub[match(nodes, rownames(fl_sub@data)),]

## Convert factor columns to characters, because this particular data is messy
indx <- sapply(fl_sub@data, is.factor)
fl_sub@data[indx] <- lapply(
    fl_sub@data[indx], function(x) as.numeric(as.character(x))
)

## Set up for compactness
fl_sub_sf<- st_as_sf(fl_sub)
pop_mat <- fl_sub_sf$pop*t(matrix(rep(fl_sub_sf$pop,nrow(fl_sub_sf)),nrow(fl_sub_sf)))
centroids <- st_geometry(st_centroid(x = fl_sub_sf))
distances <- st_distance(centroids, centroids)^2

## ---------------
## Run simulations
## ---------------
min_ssd <- 58868679
seg_truth <- read_csv("../data/largemap_enum_all/truth_fh_seg_par.csv")
seg_truth$fh <- seg_truth$fh / min_ssd
pop <- fl_sub@data$pop
rep <- fl_sub@data$mccain
nsims <- 250000
nchains <- 8
set.seed(38637) ## Random.org

## No constraints
seg_unc <- unlist(
    mclapply(1:nchains, function(x){
        out_unc <- redist.mcmc(fl_sub, popvec = pop,
                               nsims = nsims, ndists = 2)
        return(redist.segcalc(out_unc, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_unc <- data.frame(seg = seg_unc, distribution = "MCMC",
                      pop_parity = "No Equal Population Constraint",
                      compact = "No Compactness Constraint")

## Population constraints, no compactness constraint
seg_05_00 <- unlist(
    mclapply(1:nchains, function(x){
        out_05 <- redist.mcmc(fl_sub, popvec = pop,
                              nsims = nsims, ndists = 2,
                              constraint = "population", constraintweight = 10)
        out_05 <- ipw(out_05, pop = .05)
        return(redist.segcalc(out_05, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_05_00<- data.frame(seg = seg_05_00, distribution = "MCMC",
                       pop_parity = "5% Equal Population Constraint",
                       compact = "No Compactness Constraint")

seg_01_00 <- unlist(
    mclapply(1:nchains, function(x){
        out_01 <- redist.mcmc(fl_sub, popvec = pop,
                              nsims = nsims, ndists = 2,
                              constraint = "population", constraintweight = 30)
        out_01 <- ipw(out_01, pop = .01)
        return(redist.segcalc(out_01, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_01_00<- data.frame(seg = seg_01_00, distribution = "MCMC",
                       pop_parity = "1% Equal Population Constraint",
                       compact = "No Compactness Constraint")

## Compactness constraint, no population constraint
seg_00_25 <- unlist(
    mclapply(1:nchains, function(x){
        out <- redist.mcmc(fl_sub, popvec = pop,
                           nsims = nsims, ndists = 2,
                           constraint = "compact", constraintweight = .001,
                           ssd_denom = min_ssd)
        ssd_val <- unlist(lapply(1:ncol(out$partitions), function(x){
            ind <- out$partitions[, x] == 1
            return(sum(pop_mat[ind,ind]*distances[ind,ind]) +
                   sum(pop_mat[!ind,!ind]*distances[!ind,!ind]))
        })) / min_ssd
        out$compact_ssd <- ssd_val
        out <- ipw(out, com = 103827888 / min_ssd)
        return(redist.segcalc(out, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_00_25<- data.frame(seg = seg_00_25, distribution = "MCMC",
                       pop_parity = "No Equal Population Constraint",
                       compact = "25% Compactness Constraint")

seg_00_75 <- unlist(
    mclapply(1:nchains, function(x){
        out <- redist.mcmc(fl_sub, popvec = pop,
                           nsims = nsims, ndists = 2,
                           constraint = "compact", constraintweight = .01,
                           ssd_denom = min_ssd)
        ssd_val <- unlist(lapply(1:ncol(out$partitions), function(x){
            ind <- out$partitions[, x] == 1
            return(sum(pop_mat[ind,ind]*distances[ind,ind]) +
                   sum(pop_mat[!ind,!ind]*distances[!ind,!ind]))
        })) / min_ssd
        out$compact_ssd <- ssd_val
        out <- ipw(out, com = 84676202 / min_ssd)
        return(redist.segcalc(out, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_00_75<- data.frame(seg = seg_00_75, distribution = "MCMC",
                       pop_parity = "No Equal Population Constraint",
                       compact = "75% Compactness Constraint")

## 25% compact/5% population
seg_05_25 <- unlist(
    mclapply(1:nchains, function(x){
        out <- redist.mcmc(fl_sub, popvec = pop,
                           nsims = nsims, ndists = 2,
                           constraint = c("population", "compact"),
                           constraintweights = c(10, .001),
                           ssd_denom = min_ssd)
        ssd_val <- unlist(lapply(1:ncol(out$partitions), function(x){
            ind <- out$partitions[, x] == 1
            return(sum(pop_mat[ind,ind]*distances[ind,ind]) +
                   sum(pop_mat[!ind,!ind]*distances[!ind,!ind]))
        })) / min_ssd
        out$compact_ssd <- ssd_val
        out <- ipw(out, pop = .05, com = 103827888 / min_ssd)
        return(redist.segcalc(out, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_05_25<- data.frame(seg = seg_05_25, distribution = "MCMC",
                       pop_parity = "5% Equal Population Constraint",
                       compact = "25% Compactness Constraint")

## 25% compact/1% population
seg_01_25 <- unlist(
    mclapply(1:nchains, function(x){
        out <- redist.mcmc(fl_sub, popvec = pop,
                           nsims = nsims, ndists = 2,
                           constraint = c("population", "compact"),
                           constraintweights = c(30, .001),
                           ssd_denom = min_ssd)
        ssd_val <- unlist(lapply(1:ncol(out$partitions), function(x){
            ind <- out$partitions[, x] == 1
            return(sum(pop_mat[ind,ind]*distances[ind,ind]) +
                   sum(pop_mat[!ind,!ind]*distances[!ind,!ind]))
        })) / min_ssd
        out$compact_ssd <- ssd_val
        out <- ipw(out, pop = .01, com = 103827888 / min_ssd)
        return(redist.segcalc(out, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_01_25<- data.frame(seg = seg_01_25, distribution = "MCMC",
                       pop_parity = "1% Equal Population Constraint",
                       compact = "25% Compactness Constraint")

## 75% compact/5% population
seg_05_75 <- unlist(
    mclapply(1:nchains, function(x){
        out <- redist.mcmc(fl_sub, popvec = pop,
                           nsims = nsims, ndists = 2,
                           constraint = c("population", "compact"),
                           constraintweights = c(10, .01),
                           ssd_denom = min_ssd)
        ssd_val <- unlist(lapply(1:ncol(out$partitions), function(x){
            ind <- out$partitions[, x] == 1
            return(sum(pop_mat[ind,ind]*distances[ind,ind]) +
                   sum(pop_mat[!ind,!ind]*distances[!ind,!ind]))
        })) / min_ssd
        out$compact_ssd <- ssd_val
        out <- ipw(out, pop = .05, com = 84676202 / min_ssd)
        return(redist.segcalc(out, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_05_75<- data.frame(seg = seg_05_75, distribution = "MCMC",
                       pop_parity = "5% Equal Population Constraint",
                       compact = "75% Compactness Constraint")

## 75% compact/1% population
seg_01_75 <- unlist(
    mclapply(1:nchains, function(x){
        out <- redist.mcmc(fl_sub, popvec = pop,
                           nsims = nsims, ndists = 2,
                           constraint = c("population", "compact"),
                           constraintweights = c(30, .01),
                           ssd_denom = min_ssd)
        ssd_val <- unlist(lapply(1:ncol(out$partitions), function(x){
            ind <- out$partitions[, x] == 1
            return(sum(pop_mat[ind,ind]*distances[ind,ind]) +
                   sum(pop_mat[!ind,!ind]*distances[!ind,!ind]))
        })) / min_ssd
        out$compact_ssd <- ssd_val
        out <- ipw(out, pop = .01, com = 84676202 / min_ssd)
        return(redist.segcalc(out, grouppop = rep, fullpop = pop))
    }, mc.cores = detectCores())
)
seg_01_75<- data.frame(seg = seg_01_75, distribution = "MCMC",
                       pop_parity = "1% Equal Population Constraint",
                       compact = "75% Compactness Constraint")

## -----------------------------
## Create plot (without RSG yet)
## -----------------------------
truth_unc <- seg_truth %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "No Equal Population Constraint",
           compact = "No Compactness Constraint")
truth_05_00 <- seg_truth %>%
    filter(parity <= .05) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "5% Equal Population Constraint",
           compact = "No Compactness Constraint")
truth_01_00 <- seg_truth %>%
    filter(parity <= .01) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "1% Equal Population Constraint",
           compact = "No Compactness Constraint")
truth_00_25 <- seg_truth %>%
    filter(fh <= 103827888 / min_ssd) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "No Equal Population Constraint",
           compact = "25% Compactness Constraint")
truth_00_75 <- seg_truth %>%
    filter(fh <= 84676202 / min_ssd) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "No Equal Population Constraint",
           compact = "75% Compactness Constraint")
truth_05_25 <- seg_truth %>%
    filter(fh <= 103827888 / min_ssd & parity <= .05) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "5% Equal Population Constraint",
           compact = "25% Compactness Constraint")
truth_05_75 <- seg_truth %>%
    filter(fh <= 84676202 / min_ssd & parity <= .05) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "5% Equal Population Constraint",
           compact = "75% Compactness Constraint")
truth_01_25 <- seg_truth %>%
    filter(fh <= 103827888 / min_ssd & parity <= .01) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "1% Equal Population Constraint",
           compact = "25% Compactness Constraint")
truth_01_75 <- seg_truth %>%
    filter(fh <= 84676202 / min_ssd & parity <= .01) %>%
    dplyr::select(seg) %>%
    mutate(distribution = "Truth",
           pop_parity = "1% Equal Population Constraint",
           compact = "75% Compactness Constraint")

## Plot df
df_plot <- bind_rows(seg_unc,
                     seg_05_00, seg_01_00,
                     seg_00_25, seg_00_75,
                     seg_05_25, seg_05_75,
                     seg_01_25, seg_01_75,
                     truth_unc,
                     truth_05_00, truth_01_00,
                     truth_00_25, truth_00_75,
                     truth_05_25, truth_05_75,
                     truth_01_25, truth_01_75) %>%
    mutate(pop_parity = fct_relevel(pop_parity, "No Equal Population Constraint",
                                    "5% Equal Population Constraint",
                                    "1% Equal Population Constraint"),
           compact = fct_relevel(compact, "No Compactness Constraint",
                                "25% Compactness Constraint",
                                "75% Compactness Constraint"),
           distribution = fct_relevel(distribution, "Truth", "MCMC"))

## --------------
## Add RSG runs:
## --------------

# Calculate compacntess
pop <- fl_sub$pop
pop <- fl_sub$pop*t(matrix(rep(fl_sub$pop,nrow(fl_sub)),nrow(fl_sub)))
centroids <- st_geometry(st_centroid(x = fl_sub))
distances <- st_distance(centroids, centroids)^2

## UNCONSTRAINED
rsg_seg_unc <- readRDS(file = '/n/imai_lab/ckenny/enum_all/rsg_seg_unc_1.Rda')
rsg_out_unc <- readRDS(file = '/n/imai_lab/ckenny/enum_all/rsg_out_unc_1.Rda')

rsg_fh_unc <-  unlist(mclapply(1:ncol(rsg_out_unc), function(x){
    ind <- rsg_out_unc[, x] == 1
    return(sum(pop[ind,ind]*distances[ind,ind]) + sum(pop[!ind,!ind]*distances[!ind,!ind]))
}, mc.cores = detectCores() - 1))

for(ind in 2:10){
    temp_rsg_seg_unc <- readRDS(file = paste0('/n/imai_lab/ckenny/enum_all/rsg_seg_unc_', ind, '.Rda'))
    temp_rsg_out_unc <- readRDS(file = paste0('/n/imai_lab/ckenny/enum_all/rsg_out_unc_', ind, '.Rda'))
    
    temp_rsg_fh_unc <-  unlist(mclapply(1:ncol(temp_rsg_out_unc), function(x){
        ind <- temp_rsg_out_unc[, x] == 1
        return(sum(pop[ind,ind]*distances[ind,ind]) + sum(pop[!ind,!ind]*distances[!ind,!ind]))
    }, mc.cores = detectCores() - 1))
    
    rsg_seg_unc <- c(rsg_seg_unc, temp_rsg_seg_unc)
    rsg_out_unc <- cbind(rsg_out_unc, temp_rsg_out_unc)
    rsg_fh_unc <- c(rsg_fh_unc, temp_rsg_fh_unc)
    print(ind)
}

## 05
rsg_seg_05 <- readRDS(file = '/n/imai_lab/ckenny/enum_all/rsg_seg_05_1.Rda')
rsg_out_05 <- readRDS(file = '/n/imai_lab/ckenny/enum_all/rsg_out_05_1.Rda')

rsg_fh_05 <-  unlist(mclapply(1:ncol(rsg_out_05), function(x){
    ind <- rsg_out_05[, x] == 1
    return(sum(pop[ind,ind]*distances[ind,ind]) + sum(pop[!ind,!ind]*distances[!ind,!ind]))
}, mc.cores = detectCores() - 1))

for(ind in 2:10){
    temp_rsg_seg_05 <- readRDS(file = paste0('/n/imai_lab/ckenny/enum_all/rsg_seg_05_', ind, '.Rda'))
    temp_rsg_out_05 <- readRDS(file = paste0('/n/imai_lab/ckenny/enum_all/rsg_out_05_', ind, '.Rda'))
    
    temp_rsg_fh_05 <-  unlist(mclapply(1:ncol(temp_rsg_out_05), function(x){
        ind <- temp_rsg_out_05[, x] == 1
        return(sum(pop[ind,ind]*distances[ind,ind]) + sum(pop[!ind,!ind]*distances[!ind,!ind]))
    }, mc.cores = detectCores() - 1))
    
    rsg_seg_05 <- c(rsg_seg_05, temp_rsg_seg_05)
    rsg_out_05 <- cbind(rsg_out_05, temp_rsg_out_05)
    rsg_fh_05 <- c(rsg_fh_05, temp_rsg_fh_05)
    print(ind)
}

## 01
rsg_seg_01 <- readRDS(file = '/n/imai_lab/ckenny/enum_all/rsg_seg_01_1.Rda')
rsg_out_01 <- readRDS(file = '/n/imai_lab/ckenny/enum_all/rsg_out_01_1.Rda')

rsg_fh_01 <-  unlist(mclapply(1:ncol(rsg_out_01), function(x){
    ind <- rsg_out_01[, x] == 1
    return(sum(pop[ind,ind]*distances[ind,ind]) + sum(pop[!ind,!ind]*distances[!ind,!ind]))
}, mc.cores = detectCores() - 1))

for(ind in 2:10){
    temp_rsg_seg_01 <- readRDS(file = paste0('/n/imai_lab/ckenny/enum_all/rsg_seg_01_', ind, '.Rda'))
    temp_rsg_out_01 <- readRDS(file = paste0('/n/imai_lab/ckenny/enum_all/rsg_out_01_', ind, '.Rda'))
    
    temp_rsg_fh_01 <-  unlist(mclapply(1:ncol(temp_rsg_out_01), function(x){
        ind <- temp_rsg_out_01[, x] == 1
        return(sum(pop[ind,ind]*distances[ind,ind]) + sum(pop[!ind,!ind]*distances[!ind,!ind]))
    }, mc.cores = detectCores() - 1))
    
    rsg_seg_01 <- c(rsg_seg_01, temp_rsg_seg_01)
    rsg_out_01 <- cbind(rsg_out_01, temp_rsg_out_01)
    rsg_fh_01 <- c(rsg_fh_01, temp_rsg_fh_01)
    print(ind)
}

# Put everything together for rsg
rsg_df_unc <- tibble(seg = rsg_seg_unc, distribution = 'RSG', 
                     pop_parity = "No Equal Population Constraint",
                     fh = rsg_fh_unc)

rsg_df_05 <- tibble(seg = rsg_seg_05, distribution = 'RSG', 
                    pop_parity = "5% Equal Population Constraint",
                    fh = rsg_fh_05)

rsg_df_01 <- tibble(seg = rsg_seg_01, distribution = 'RSG', 
                    pop_parity = "1% Equal Population Constraint",
                    fh = rsg_fh_01)


df_plot_rsg <- bind_rows(rsg_df_unc, rsg_df_05, rsg_df_01)

# create subsets to match compactness constraints
df_plot_rsg <- bind_rows(
    (df_plot_rsg %>% filter(fh <= 84676202) %>% mutate(compact = '75% Compactness Constraint')),
    (df_plot_rsg %>% filter(fh <= 103827888) %>% mutate(compact = '25% Compactness Constraint')),
    (df_plot_rsg %>% mutate(compact = 'No Compactness Constraint'))
) %>% select(-fh)

## Join run types and plot:
df_plot <- bind_rows(df_plot, df_plot_rsg)

df_plot <- df_plot %>%
    mutate(pop_parity = fct_relevel(pop_parity, "No Equal Population Constraint",
                                    "5% Equal Population Constraint",
                                    "1% Equal Population Constraint"),
           compact = fct_relevel(compact, "No Compactness Constraint",
                                 "25% Compactness Constraint",
                                 "75% Compactness Constraint"),
           distribution = fct_relevel(distribution, "Truth", "MCMC", "RSG"))

df_plot <- df_plot %>% mutate(compact = fct_recode(compact, "25th Percentile\n Compactness Constraint" = "25% Compactness Constraint",
                            ),
                           compact = fct_recode(compact, "75th Percentile\n Compactness Constraint" = "75% Compactness Constraint",
                           ))


ggplot(df_plot, aes(seg, colour = distribution, fill = distribution,
                    alpha = distribution, lty = distribution)) +
    geom_density(lwd = 1.1) +
    facet_grid(compact ~ pop_parity) +
    theme_classic() +
    scale_colour_manual(values = c("grey", "black", "red")) +
    scale_fill_manual(values = c("grey", "white", "white")) +
    scale_linetype_manual(values = c(1, 1, 2)) + 
    scale_alpha_manual(values = c(1, 0, 0)) + 
    labs(x = "Republican Dissimilarity Index", y = "Density") +
    ben_theme() +
    ggsave("largemap_enum_all_compact_validation.pdf", height = 8, width = 8)
    #ggsave("../paper/figs/largemap_enum_all_validation.pdf", height = 10, width = 10)

## ## --------
## ## RSG sims
## ## --------

## ## Load RSG sims
## load("../data/largemap_enum_all/simulations_rsg.RData")

## rsg_unc <- do.call("cbind", lapply(sims_out, "[[", 1))
## seg_unc_rsg <- redist.segcalc(rsg_unc, grouppop = rep, fullpop = pop)

## rsg_05 <- do.call("cbind", lapply(sims_out, "[[", 2))
## seg_05_rsg <- redist.segcalc(rsg_05, grouppop = rep, fullpop = pop)

## rsg_01 <- do.call("cbind", lapply(sims_out, "[[", 3))
## seg_01_rsg <- redist.segcalc(rsg_01, grouppop = rep, fullpop = pop)

## -------------------------------
## Descriptive plot of the new map
## -------------------------------

## Reset some objects to how they're initially loaded for easier use
seg_truth <- read_csv("/n/imai_lab/ckenny/enum_all/truth_fh_seg_par.csv")

fl_map <- readShapePoly("../data/fl/FL.shp")

## Load nodes
mn <- 89 
nodes <- strsplit(readLines("../data/largemap_enum_all/map_sub_nodes_89.dat")," ")[[1]]

## First, subset down to right map
fl_sub <- fl_map[which(rownames(fl_map@data) %in% nodes),]

## Then, reorder nodes to be correct order to match - otherwise, won't
## line up properly with original map
fl_sub <- fl_sub[match(nodes, rownames(fl_sub@data)),]

## Convert factor columns to characters, because this particular data is messy
indx <- sapply(fl_sub@data, is.factor)
fl_sub@data[indx] <- lapply(
    fl_sub@data[indx], function(x) as.numeric(as.character(x))
)


## Begin analysis
fl_sub@data$id <- rownames(fl_sub@data)
fl_sub_points <- fortify(fl_sub, region = "id")
fl_sub_df <- plyr::join(fl_sub_points, fl_sub@data, by = "id")

## Map plot
map_plot <- ggplot(fl_sub_df) + 
    aes(long, lat, group = group, fill = pop) + 
    geom_polygon(colour = "black", size = .3) +
    geom_path(color="black", lwd = .05) + 
    theme_classic() +
    scale_fill_gradient(low = "white", high = "black") +
    coord_equal() +
    guides(fill = guide_legend(title = "Population")) +
    labs(x = "Longitude", y = "Latitude",
         title = "70-Precinct Validation Map") + 
    theme(panel.border = element_rect(colour = "black", fill=NA),
          plot.title = element_text(hjust = 0.5),
          legend.position = c(.11, .15),
          legend.title = element_text(size = 8), 
          legend.text = element_text(size = 8))

## Parity plot
### No compactness constraint
parity_plot_df <- seg_truth %>% filter(parity <= .2) %>%
    mutate(
        parity_group = cut(parity, breaks = seq(0, max(parity), .01),
                           include.lowest= TRUE, labels = FALSE)
    ) %>%
    group_by(pg = parity_group / 100) %>%
    summarize(n = n())
numsols_df <- c(
    sum(seg_truth$parity <= .01),
    sum(seg_truth$parity <= .05),
    sum(seg_truth$parity <= .1)
)
text_out <- data.frame(text = paste0(paste0("< ", format(c(.01, .05, .10), nsmall = 2), " Parity: ", numsols_df), collapse = "\n"))

parity_plot <- ggplot(parity_plot_df, aes(x = pg - .005, y = n)) +
    geom_bar(stat = "identity", orientation = 'x') +
    theme_classic() +
    scale_y_continuous(labels = scales::comma) +
    geom_text(data = text_out, aes(x = .1, y = Inf, label = text),
              vjust = 1.2, hjust = 1.1) + 
    labs(x = "Distance from Population Parity",
         y = "Number of Plans",
         title = "Distribution of Population Parity on 70-Precinct Map\n with No Compactness Constraint") + 
    theme(panel.border = element_rect(colour = "black", fill=NA),
          plot.title = element_text(hjust = 0.5),
          legend.position = c(.05, .05)) + 
    xlim(0, .20)+
    ylim(0, 1200000)


### 25 compactness constraint
parity_plot_df_25 <- seg_truth %>% filter(parity <= .2, fh <= 103827888) %>%
    mutate(
        parity_group = cut(parity, breaks = seq(0, max(parity), .01),
                           include.lowest= TRUE, labels = FALSE)
    ) %>%
    group_by(pg = parity_group / 100) %>%
    summarize(n = n())
numsols_df_25 <- c(
    sum(seg_truth$parity <= .01 & seg_truth$fh <= 103827888),
    sum(seg_truth$parity <= .05 & seg_truth$fh <= 103827888),
    sum(seg_truth$parity <= .1 & seg_truth$fh <= 103827888)
)
text_out_25 <- data.frame(text = paste0(paste0("< ", format(c(.01, .05, .10), nsmall = 2), 
                                            " Parity: ", numsols_df_25), 
                                     collapse = "\n"))

parity_plot_25 <- ggplot(parity_plot_df_25, aes(pg - .005, n)) +
    geom_bar(stat = "identity", orientation = 'x') +
    theme_classic() +
    scale_y_continuous(labels = scales::comma) +
    geom_text(data = text_out_25, aes(x = .1, y = Inf, label = text),
              vjust = 1.2, hjust = 1.1) + 
    labs(x = "Distance from Population Parity",
         y = "Number of Plans",
         title = "Distribution of Population Parity on 70-Precinct Map\n with 25th Percentile Compactness Constraint") + 
    theme(panel.border = element_rect(colour = "black", fill=NA),
          plot.title = element_text(hjust = 0.5),
          legend.position = c(.05, .05))+
    xlim(0, .20) + 
    ylim(0, 1200000)


### 75 copmactness constraint
parity_plot_df_75 <- seg_truth %>% filter(parity <= .2, fh <= 84676202) %>%
    mutate(
        parity_group = cut(parity, breaks = seq(0, max(parity), .01),
                           include.lowest= TRUE, labels = FALSE)
    ) %>%
    group_by(pg = parity_group / 100) %>%
    summarize(n = n())
numsols_df_75 <- c(
    sum(seg_truth$parity <= .01 & seg_truth$fh <= 84676202),
    sum(seg_truth$parity <= .05 & seg_truth$fh <= 84676202),
    sum(seg_truth$parity <= .1 & seg_truth$fh <= 84676202)
)
text_out_75 <- data.frame(text = paste0(paste0("< ", format(c(.01, .05, .10), nsmall = 2), 
                                               " Parity: ", numsols_df_75), collapse = "\n"))

parity_plot_75 <- ggplot(parity_plot_df_75, aes(pg - .005, n)) +
    geom_bar(stat = "identity", orientation = 'x') +
    theme_classic() +
    scale_y_continuous(labels = scales::comma) +
    geom_text(data = text_out_75, aes(x = .1, y = Inf, label = text),
              vjust = 1.2, hjust = 1.1) + 
    labs(x = "Distance from Population Parity",
         y = "Number of Plans",
         title = "Distribution of Population Parity on 70-Precinct Map\n with 75th Percentile Compactness Constraint") + 
    theme(panel.border = element_rect(colour = "black", fill=NA),
          plot.title = element_text(hjust = 0.5),
          legend.position = c(.05, .05)) +
    xlim(0, .20) + 
    ylim(0, 1200000)


pdf(file = "../paper/figs/map_descriptives_70precinct.pdf", height = 10, width = 10)
grid.arrange(map_plot, parity_plot, parity_plot_25, parity_plot_75, widths = c(2, 2), nrow = 2)
dev.off()

library(patchwork)
(map_plot + parity_plot)/(parity_plot_25 + parity_plot_75)
ggsave("../paper/figs/map_descriptives_70precinct.png")
